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We use large deviation methods to calculate rates of noise-induced transitions between states in 
multi-stable genetic networks. We analyze a synthetic biochemical circuit, the toggle switch, and 
compare the results to those obtained from a numerical solution of the master equation. 
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Fluctuations in bio-molecular networks have been the 
subject of much research activity recently ll|. Studies 
on noise in gene expression J2|, y, U, la, la, |j], i n signal 
transduction |a| and in biochemical oscillators |3, |ljj, |jjj 

demonstrated that having a small number of molecules 
affects, sometimes critically, the behavior of cellular cir- 
cuits. Stochastic aspects of the choice between lytic 
and lysogenic developmental strategies of bacteriophage 
lambda virus infection in E. coli were studied in an in- 
fluential paper by Arkin, Ross and McAdams |l2(. 

One of the interesting aspects developmental processes 
is that one could get multiple heritable cell fates without 
irreversible changes to the genetic information. Differ- 
ent cells with the same DNA sequence, showing different 
phenotypes that are stably maintained through cell di- 
visions, namely epigenetic phenomena, have been rep- 
resented as multiple stable attractors in deterministic 
descriptions of the biochemical dynamics. In this pa- 
per, we are concerned with the robustness of such at- 
tractors against spontaneous fluctuations which might 
induce transitions from one stable state to another. Pre- 
vious work in this area has modeled the effects of fluc- 
tuations by adding Gaussian-distributed Langevin forces 
to the deterministic equations [iji LLJ, LL£| . Although this 
description is appropriate in describing typical fluctua- 
tions when the number of molecules is sufficiently large 
|2L VA 111 la- Tare events involving occasional large de- 
partures from average behavior are typically outside the 
scope of the Langevin treatment (Gaussian approxima- 
tion). The transition rate in a simplified model of the 
phage lambda switch has been studied [1J, |lj| in this 
approximation. We wish to compute the transition rate 
using a more appropiate large deviation theory with spe- 
cial focus on the attempt frequency. Of course, one could 
get the transition rate from direct computer simulations. 
However, direct simulations of rare events is, obviously, 



time-consuming. Recent research in the lambda switch 
suggests that the simplified model lacks one very im- 
portant physical interaction between distant regions of 
the lambda virus genome, changing dramatically the be- 
haviour of the switch |lj|. Applying our tools to that 
question, among others, is the long time goal of our re- 
search. However, we wish to test our methods on a sim- 
pler system. We will consider the artificially constructed 
toggle switch |l8|. In this example, we find that the con- 
tributions to the transition rate coming from corrections 
to the Gaussian approximation can change the overall 
rate by several orders of magnitude and, therefore, are 
important for comparison with experimental results. 

The theory of transition rates is a well developed sub- 
ject (see jlSj as well as references therein). For a bistable 
system like the genetic switch we are considering, the 
transition probability from one stable point to the other 
is estimated by computing the probability of reaching 
the saddle point between stable states, and, from there, 
to follow the deterministic trajectory to the other stable 
state, rather than to fall back to the initial state. The 
transition rate is given by an expression of the form: 
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where A+ is the positive eigenvalue of the matrix describ- 
ing the linearized equations of motion around the saddle 
point, Af p and A sp are the inverses of covariance matrices 
appearing in the quasi-stationary Gaussian approxima- 
tion of the probability distribution in the starting stable 
point and in the saddle point respectively and P(xf,x ) 
is the probability of finding the system at the saddle point 
state in a quasi-stationary distribution centered around 
the stable fixed point x . Note that A sp has one nega- 
tive eigenvalue and the Gaussian distribution around the 
saddle point is only a formal solution. A more precise 



definition of A appears later in the pap er. For a deriva- 
tions of a very similar formula see |20| or section VII. D 
in the review |l9j . 

Much of the rest of the paper is devoted to the com- 
putation of P{xf 1 x ) by large deviation methods. There 
are two related ways. In one approach, one keeps track 
of the trajectories in the space of numbers of different 
molecules, distributed according to a state dependent 
Poisson process, and computes time dependent transi- 
tion probabilities as sum the probabilities of all paths 
connecting the initial and final points, which leads nat- 
urally to a path integral formulation of the stochastic 
process. In this way, the transition probability is eval- 
uated as the exponential of the "action". This action 
can be computed in a perturbation expansion (using the 
volume of the system as a parameter) , in which the lead- 
ing order correction is the line integral along the path 
that minimizes the action (optimal path) of a Lagrangian 
function. This calculation naturally gives rise to a Hamil- 
tonian that corresponds to the evolution operator in the 
master equation [2l|, written in terms of numbers and 
raising and lowering operators expressed as exponentials 
of the phase variables conjugate to the numbers. 

An alternative but exactly equivalent approach is to 
start directly from the master equation and solve it in the 
Eikonal approximation 22]. We will present our argu- 
ments in this article using this approach, which is easier 
to explain mathematically, and provides an easier way to 
compute the next order correction in the volume expan- 
sion, a term which has not been computed before in rela- 
tion to these genetic switches. As we will see, in the case 
of the toggle switch, the next order term in In P{xf,x ) 
makes an important contribution to the overall rate of 
transition. 

The general ideas are developed in the context of the 
simple example of the toggle switch. This artificially re- 
alized switch consists of two genes that repress each oth- 
ers' expression, placed in a high copy plasmid in E. coli. 
Once expressed, each protein can bind particular DNA 
sites upstream of the gene which codes for the other pro- 
tein, thereby repressing its transcription. If we denote 
the i-th protein concentration by Xi, the deterministic 
system is described by the equations: 

a% x\ 
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to and n, represent the degree of cooperative binding 
of proteins to DNA, and t _1 is the protein degrada- 
tion/dilution rate (assumed equal for the two proteins). 
K\ is the effective dissociation constant for binding of 
protein 1 in the promoter of gene 2. K% is the corre- 
sponding parameter for protein 2. For some regions of 
parameter space, the system has three stationary points: 
two stable ones and a saddle point |l8| ■ 

For the purposes of this discussion, we model the 
stochastic evolution of the protein concentrations in the 
system by a birth-death process in which protein i is 
made in short-lived bursts of size bi and proteins are di- 
luted or degraded at a rate r — 1 . A more detailed de- 
scription involving proteins and RNA will be published 
elsewhere. It is worth noting that, while both the burst 
size bi and the RNA production rate show up as param- 
eters in the stochastic modeling, only their product, Oj, 
shows up in the effective deterministic equations (J2J) for 
the protein levels. 

To compute the rate of transition from one fixed point 
to the other, we must solve the master equation [21|, 
which describes the time evolution of the probability dis- 
tribution of protein concentrations. The qualitative be- 
havior of the stationary solution for the bistable system 
can be described in simple intuitive terms: the solution 
displays two peaks centered around the stable points. If 
we start with probability one around one of the stable 
points, rare transitions lead to a long tail which leaks 
into the domain of attraction of the other stable point, 
in very much the same way in which the probability am- 
plitude extends beyond the classically allowed region in 
quantum mechanical tunneling through a barrier. This 
analogy motivates the Eikonal approximation to the so- 
lution of the master equation [22| . The master equation 
is given by, 

OP 

— = nY^[Wg{x - i/n)p{x - e/n,t) - Wi(x)P(x,t)] 

e 

(4) 
where Q is the volume of the system, e/f2 = Ax is the 
concentration change associated with individual reaction 
events, the rate of which is given by QWe(x). Assuming 
that the distribution is quasi-stationary in the region of 
interest, we consider solutions of the WKB form: 



P(x,t)=Cexp[-nS(x),] S(x o ) = 0. 



(5) 



1 + {xx/KxY 

the constants a\ and a-i incorporate all aspects of tran- 
scription and translation reactions. The Hill exponents 



x being the initial stable point. In the same way the 
wave function in quantum mechanics is computed using 
an expansion in powers of h, it customary to find the 



probability P(x,t) by expanding S(x) in powers of in- 
verse volume, which plays the same role as h in quantum 
mechanics, since the bigger the volume, the less likely are 
fluctuations to happen. Then, to first order in fi _1 , we 
write: 

S(x) = S Q (x) + fi-^ifas) + 0(£T 2 ). 

Assuming that the scaled transition rates W&(x) are 
smooth functions of x, and expanding S to first order, 
S(x — e/Q) = S(x) — j£.-7^-S(x), collecting the terms 
which do not contain powers of $7 we have: 



dP(x,t) 
di 



HP(x,t) 



H(x,p)=Y / [We(x)(e^^-l) 



(6) 

(7) 



used in [13|, |l4|. Since we are considering a situation 
where the transitions are so rare that the probability does 
not change much in time, the Hamiltonian will be very 
small. 

The main contribution to the transition probability is 
obtained by evaluating P along a particular trajectory 
|22j. This trajectory, called the optimal path, is the so- 
lution to Hamilton's equations derived from Eq0 
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where H is the Hamiltonian describing the time evolu- 
tion of the probability distribution, and we define the 
momentum pi as: 



Pi = t^—S q (x) 



(8) 



If we expand the Hamiltonian0in p and keep terms up 
to second order in p we recover the Gaussian approach 



For the toggle switch example we have four e^-s de- 
scribing jumps to the right, left, up or down, given by 
bxXx, — xx,b%Xi, and — £2, respectively. The relevant 
Hamiltonian defined on times long compared to the in- 
verse binding/unbinding rates of proteins at the two pro- 
moters is given by: 
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As already mentioned above, Kx,i are the effective disso- 
ciation constants for binding of proteins 1,2 at the pro- 
moter of gene 2, 1, respectively, bi is the burst size of 
protein i and the ratio ai/bi is a measure of the RNA 
production rate associated with the transcription of the 
gene i. 

To extract the values of the burst size parameters, the 
spontaneous transition rate has to measured experimen- 
tally for more than one conditions. Since this has not 
yet been done, we will compare the results of the Eikonal 
approximation to the solution obtained by direct diago- 
nalization of the Hamiltonian (jl 1|> . For simplicity we will 
set the parameters Ki = 1,6, = 1. 

The optimal path for the transition from one stable 
point to the other starts near one stable point, proceeds 
to the saddle point and from there it follows the deter- 
ministic trajectory to the other stable point. Thus we 



must first find solutions of Eqs. I§1 and ITUl which start at 
(near) the initial stable point and end at the saddle point. 
At the end points we have px = P2 = 0, and H = 0. This 
also implies that if the system is at the stable point it 
will remain there. So, the optimal path must instead 
start at a point very close to but not exactly at the fixed 
point. In this case, the Hamiltonian will be a very small 
number (and constant). In what follows, we will make 
the approximation H — 0. The initial conditions for the 
momentum equations can be obtained by approximating 
the probability around the stable point by a Gaussian dis- 
tribution P = e~ nSg with S g = ^AijSxiSxj (note that 
we use summation convention, i.e., repeated indices are 
summed over). Then pi = ^rr = AijSxj, and we expand 
the equation H = around the stable point to find Aij . 
Then we have a two point boundary value problem which 
can be solved by various methods [2j|. The solution of 



the equations of motion I§1 and ITU1 for a set of parameters, 
projected to concentration space, is shown in Figure Q 
We integrate equations |H1 along the optimal path C to 
obtain Sq = J c Pidxi . 
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FIG. 1: Optimal path for the parameters, ai = 156, a,2 — 
30,n = 3,m = 1, K\ — K2 — 1, 61 = 62 = 1 and r = 1. xi 
are dimensionless. The ellipsoid indicates the orientation of 
the Gaussian spread around the stable point. The size of the 
spread scales like fi~2. 

The S\ factor can be viewed as a correction due to 
fluctuations around the optimal path and could be cal- 
culated following references J2J, l25l |. Collecting coeffi- 
cients of powers of Q, in the fl _1 expansion we derive an 
equation for Si: 



^2[Wsei^- - ^eiijdiPj ~ eAWi\e^ Va) = (12) 

e 

In turn, after using the equations of motion to rewrite 
the first term as derivative along the optimal path x op (t') 1 
Eq. (|12l) can be transformed into: 
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(13) 



To proceed we need j^f- along the path. From Hamil- 
ton's equations © it follows that 5p(t) a = M(t) a bSx(t)b, 
and thus we can use the components of the matrix M in 
place of the derivative -g£ in (|12fl . Moreover, J5J also 
implies that: 

d 2 H . i d 2 H . 
Sx = a » i Sx + a—a- S P* 14 

.. d 2 H . d 2 H 

SPa = —7; 7T~" X - "?; TT~ "Pi (15) 

ax a ox % ax a opi 



Combining this together with the time derivative of Sp(t) , 

6p = M5x + MSx (16) 

leads to the following set of coupled differential equations 
for M: 

r) 2 H B 2 H 

M ab + M ac - j- +M aC n - M db 



'dx b dp c dp c dpd 

d 2 H w d 2 H 

-M cb + ^— .^ = 



dx a dp c 



dx a dx h 



(17) 



with initial conditions: A/y (t = 0) = Ay (defined below 
equation llljl . Finally, solving these equations together 
with equations |H| and ^J we integrate equation (|13fl to 
obtain Si . Given the above values of So and Si we com- 
pute the transition probability, P{xf, x ), from the start- 
ing stable point, x , to the saddle point, Xf. Using equa- 
tion ^we can, therefore, find the transition rate for any 
large value of fL We now compare this calculation to the 
direct estimation of transition rates as described below. 
From the master equation |J3J| , it follows that the eigen- 
values of H measure the decay rates of non-stationary 
states corresponding to eigenvectors of H with nonzero 
eigenvalues. The equilibrium state is represented by the 
"zero mode", i.e., the eigenvector of H with zero eigen- 
value, the existence of which is guaranteed by the tran- 
sition matrix character of the Hamiltonian and conser- 
vation of probability. To compute the eigenvalues of the 
Hamiltonian, we write the master equation in discrete 
form, replacing the continuous concentration variables 
(xi,X2) with a lattice with lattice parameter I/O. Al- 
though the system displays infinitely many states, typi- 
cally, the gap between the real parts of the eigenvalues 
for first and second excited states is much larger than the 
absolute value of the real part of the first eigenvalue. This 
is because the gap between the first excited state and the 
second or the third excited states are governed by local 
relaxation rate around the two fixed points, but, the gap 
between the ground state and the first excited state is 
governed by the transition rate between the two stable 
fixed points. The local relaxation rates are order one in 
O, whereas, the transition rate is exponentially small for 
large f2 (in practice, we find the ratios of the real parts 
to be about 10 3 ). Thus an arbitrary probability distri- 
bution rapidly decays into a linear combination of the 
stationary state and the first excited state. Equivalently, 
the state could be described as a linear combination of 
two states, each representing a quasi-stationary distribu- 
tion around a stable fixed point. From then on, we can 



project the evolution to this two state system. If we start 
with probability p of being in the state (1, 0) T , then the 
Master equation gives: 




The two-by-two effective transition matrix has columns 
which sum to zero ensuring probability conservation. 
Also, the trace + e\ = r i2 + r 2 i, where ei is the eigen- 
value of the first excited state. Therefore the first excited 
eigenvalue will be the sum of the forward and backward 
rates. In the case of the asymmetric systems, one rate 
is usually far greater than the other. Consequently the 
larger rate among r\% and r%\ will be approximately given 
by ei, which we computed numerically using the Matlab 
routine "ergs" for sparse matrices as well as by Lanczos 
algorithm 26] . For a symmetric choice of parameters for 
the two proteins, each rate is just ei/2. 

To explicitly extract the Sq and Si contributions to the 
rate from the Lanczos results, we re-scale the volume of 
the system Q — > vQ which, in turn, leads to a re-scaling 
of rates of individual reaction events as f(x) — > vf{x). 
As a function of volume scale factor, v, the logarithm of 
the rate has the form: ln(r) = Sqis + 6, where b includes 
both Si and the logarithm of the pre- factor of P(xf,x ) 
in EqQ] The results and comparison with the Eikonal ap- 
proximation are shown in Fig|2 The dotted line is a fit 
to the data points obtained from calculation of the eigen- 
values, and we see that the slope and intercept computed 
from eauations lll5l are in good agreement with these val- 
ues. Note that, in this example, Si and the pre-factor 
are significant contributions to the transition rate. 

When we perform these calculations for the "standard" 
model of the lambda switch [14J |27j , we find a rate three 
orders of magnitude higher than the observed rate of 10 -7 
per generation |28j . In retrospect, it is clear that account- 
ing for the stability of the lysogenic state requires a more 
complex model which should include the effect of DNA 
looping [F7J. Whether the stability is due to suppression 
of fluctuation or due to disappearance of the lytic "fixed" 
point 29] remains an open question. 



Optimal path methods are routinely used for study- 
ing rare events related to failure of communication net- 
works modeled as birth and death processes [3(|. Such 
large deviation methods are likely to be important in the 
context of robustness and adaptability of biological net- 
works. This paper illustrates the power of an approach to 
fluctuations based on the Eikonal approximation to so- 
lutions of the master equation. The scheme incorporates 
-6 
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FIG. 2: Scaling with volume: estimates from direct compu- 
tation of eigenvalues are So — 2.63, Si + ln(pref) = 4.85 
whereas optimal path calculation gives So = 2.47, Si = 3.5, 
ln(pref) = 1.5. In this example the backward rate is 1000 
times smaller than the forward rate, so the lowest nonzero 
eigenvalue is very close to the rate of switching. 



large deviations in a natural way and provides a quanti- 
tative method scalable to large networks. We also hope 
that beyond being an efficient computational tool, this 
method will provide further insight into to the stability 
of epigenetic states of complex genetic networks. 
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